# S. aureus infect mice skin
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

key_cytokine <- c('Ccl20','Tslp','Flt3l','Csf2','Tnf','Il6')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg')

# realign raw data ------
qzeng.path <- list.files('/home/supervisor/mist2/gjsx/qzeng2024il24', 'gz', full.names = T)

tibble(fq = qzeng.path) |>
  mutate(sample = ifelse(str_detect(fq, '097'), 'Infected', 'PBS'),
         readgroup = sample,
         end = ifelse(str_detect(fq, '_1'), 'fastq_1', 'fastq_2')) |>
  pivot_wider(names_from = end, values_from = fq) |>
  write_csv('/home/supervisor/mist2/gjsx/qzeng2024il24/nf.starsolo.csv')

# 24311 gene 20485 cells
sobj <-
c('/home/supervisor/mist2/gjsx/qzeng2024il24/myindex_unqiue_em_out/PBS/Solo.out/GeneFull_Ex50pAS/filtered/',
  '/home/supervisor/mist2/gjsx/qzeng2024il24/myindex_unqiue_em_out/Infected/Solo.out/GeneFull_Ex50pAS/filtered/') |>
  set_names(c('PBS','infected')) |>
  Read10X(strip.suffix = TRUE) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

# matrix read in -----------
cnt <- c('mission/FPP/zww_sa_mice/data/zww-matrix-a/', 'mission/FPP/zww_sa_mice/data/zww-matrix-b/') |>
  set_names(c('PBS','infected')) |>
  Read10X(strip.suffix = TRUE)

sobj <- cnt |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj$mito.ratio <- sobj |> PercentageFeatureSet('^mt-')

sobj |> VlnPlot('mito.ratio')

sobj <- sobj |> filter(mito.ratio < 5)

sobj <- sobj |> quick_process_seurat()

# auto annotate ------
micerna <- celldex::MouseRNAseqData()

sobj <- sobj |>
  mark_cell_type_singler(ref = micerna, new_label = 'bulk_main')

DimPlot(sobj, group.by = 'bulk_main', cols = 'Paired') +
  ggtitle('Main cell types')

DimPlot(sobj, group.by = 'bulk_main', cols = my36colors, split.by = 'orig.ident') +
  ggtitle('Main cell types')

FeaturePlot(sobj, c('Krt14','Krt15'))

# get all marker ----------
all.marker <- FindAllMarkers(sobj, logfc.threshold = 1, only.pos = T)

krt_fmly <- filter(all.marker, str_detect(gene, 'Krt')) |>
  pull(gene) |>
  unique() |>
  sort()

DotPlot(sobj, c('Krt14','Krt15'))

# kera marker from panglaodb
pangl.kc <- c('Krt10','Krt15','Krt16','Ivl','Dmkn','S100a14')

sobj |> DotPlot(pangl.kc, cluster.idents = T)

sobj |> DotPlot(krt_fmly) +
  RotatedAxis()

sobj |>
  filter(bulk_main == 'Fibroblasts') |>
  DotPlot(c('Krt14','Krt15'))

sobj |>
  filter(bulk_main == 'Fibroblasts') |>
  DotPlot(krt_fmly) +
  RotatedAxis()

sobj <- sobj |>
  mutate(manual_main = case_when(
    seurat_clusters %in% c(2,20,21,6,4,3,19,18,7,5,9) ~ 'Keratinocytes',
    seurat_clusters == 11 ~ 'DC',
    .default = bulk_main
  ))

sobj |> write_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

sobj <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

### FIG: mice skin cell type umap =============
sobj |>
  DimPlot(group.by = 'manual_main', cols = DiscretePalette(36),
          label = T, label.size = 1.5) +
  ggtitle('Main cell types in mice skin') +
  theme_jpub +
  NoLegend() 

publish_pdf('micefig/mice_skin_celltype_umap.pdf')

g1 <- last_plot()

HoverLocator(g1, information = FetchData(sobj, vars = 'seurat_clusters'))

# fraction of main types ------
### FIG: DC frac change following SA infection =========
dc_conf <- sobj |>
  mutate(DC = ifelse(manual_main %in% c('Granulocytes','DC'), 'DC', 'nonDC')) |>
  dplyr::count(orig.ident, DC) |>
  calc_frac_conf_on_grouped_count()

sobj |>
  mutate(DC = ifelse(manual_main %in% c('Granulocytes','DC'), 'DC', 'nonDC')) |>
  pluck('meta.data') |>
  test_on_grouped_count(orig.ident, DC) 

dc_conf |>
  filter(DC == 'DC') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(orig.ident, fraction, ymax = conf.high, ymin = conf.low, fill = orig.ident)) +
  geom_col() +
  geom_errorbar(width = .5) +
  labs(x = 'Group', y = 'Fraction in all skin cells', fill = 'Group') +
  scale_fill_manual(values = c('skyblue','orange')) +
  theme_jpub

publish_pdf('micefig/mice_DC_frac_PBS-SA.pdf')

dc_deg <- sobj |>
  filter(manual_main %in% c('Granulocytes','DC')) |>
  FindMarkers(group.by = 'orig.ident',
              ident.1 = 'infected',
              logfc.threshold = 1,
              only.pos = T)

dc_upgo <- dc_deg |>
  filter(p_val_adj < .05) |>
  rownames() |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T) |>
  as_tibble()

dc_upgo[1:10,] |>
  publish_enrichment()

publish_pdf('micefig/mice_SA_DC_upgo.pdf', width = 60)

# focus on kera ------
sobj_kera <- sobj |>
  filter(manual_main == 'Keratinocytes')

## compare infect kera vs PBS kera -----
kera.infect.deg <- sobj_kera |>
  FindMarkers(ident.1 = 'infected', group.by = 'orig.ident',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene')

kera.infect.deg |> write_csv('mission/FPP/zww_sa_mice/results/kera.infect.deg.csv')

kera.infect.deg |>
  filter(gene %in% key_cytokine) |>
  ggplot(aes(x = 'Keratinocyte', gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  theme_jpub

sobj_kera <- sobj_kera |>
  quick_process_seurat(skip_norm = T)

sobj_kera <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

sobj_kera |> write_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

sobj_kera <- sobj_kera |>
  FindClusters(resolution = .7, algorithm = 4,
               method = 'igraph', cluster.name = 'leiden0.7')

sobj_kera |> DimPlot(cols = DiscretePalette(36), label = T)

## remove cluster with high mito
sobj_kera |>
  VlnPlot('mito.ratio', pt.size = 0, cols = DiscretePalette(36)) +
  stat_mean()
  
sobj_kera <- sobj_kera |> filter(leiden0.7 != 9)

sobj_kera$seurat_clusters <- sobj_kera$leiden0.7 |>
  fct_drop() |>
  as.numeric() |>
  as.character() |>
  fct_infreq() 

## FIG: 11 clusters kera umap =============
Idents(sobj_kera) <- 'seurat_clusters'

sobj_kera |> DimPlot(cols = DiscretePalette(36), label = T, label.size = 2) +
  theme_jpub +
  NoLegend()

publish_pdf('micefig2/mice_kera_11c_umap.pdf')

## find Trpv3-high keratinocytes ----------
### Trpv3 with late/early KC marker ----------
sobj_kera |>
  DotPlot(c('Trpv3', late.kc, 'Krt5','Krt14'), group.by = 'seurat_clusters',
          cluster.idents = T, cols = 'RdYlBu', dot.scale = 3) +
  labs(title = 'Trpv3 and KC differentiation markers in mice skin KC',
       x = 'Genes', y = 'KC clusters') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_kera_late_marker_v3.pdf', width = 55)

### FIG: PBS-Trpv3 SA-Ccl20 dotplot ============
sobj_kera |>
  filter(orig.ident == 'PBS') |>
  DotPlot(c('Trpv3'), group.by = 'seurat_clusters', dot.scale = 3) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  labs(x = 'Gene', y = 'Keratinocyte clusters', title = 'PBS group')

publish_pdf('micefig2/mice_kera_pbs-v3_bubble.pdf', width = 40)

g1 <- last_plot()

sobj_kera |>
  filter(orig.ident == 'infected') |>
  DotPlot(c('Ccl20'), group.by = 'seurat_clusters') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  labs(x = 'Gene', y = 'Keratinocyte clusters', title = 'SA group')

g2 <- last_plot()

g1 + NoLegend() + g2 
  
publish_pdf('micefig2/mice_kera_pbs-v3_SA-ccl20_bubble.pdf', width = 60)

### find Trpv3-high kera in PBS ---------
Idents(sobj_kera) <- 'seurat_clusters'

trpv3_rich_pbs <- sobj_kera |>
  filter(orig.ident == 'PBS') |>
  FindAllMarkers(features = 'Trpv3', logfc.threshold = 0)

trpv3_rich_pbs |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(cluster, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 1, linetype = 'dashed') +
  scale_fill_manual(values = c('blue','grey','red')) +
  labs(title = 'Enrichment of Trpv3 expression in uninfected keratinocytes',
       y = 'Enrichment log2FC',
       caption = 'clusters with <1% Trpv3 expression are not shown') +
  theme_pubr() + labs_pubr()

g1 <- last_plot()

### find ccl20-high kera in SA -----------
ccl20_rich_sa <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers(features = 'Ccl20', logfc.threshold = 0)

ccl20_rich_sa |>
  bind_rows(trpv3_rich) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(cluster, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 1, linetype = 'dashed') +
  scale_fill_manual(values = c('blue','grey','red')) +
  facet_wrap(~gene, ncol = 1, scales = 'free_y') +
  labs(y = 'Enrichment log2FC') +
  theme_jpub

publish_pdf('micefig2/mice_trpv3_ccl20_enrich_in_kera_barplot.pdf')

### FIG: kera PBS-trpv3 SA-ccl20 enrich bubble ===========
ccl20_rich_sa |>
  bind_rows(trpv3_rich) |>
  ggplot(aes(y = cluster, x = gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  labs(caption = 'clusters with <1% gene expression are not shown') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(breaks = c(0,75,150)) +
  theme_jpub

publish_pdf('micefig2/mice_trpv3_ccl20_enrich_in_kera_bubble.pdf')

### FIG: co-express of Ccl20 & Trpv3 in kera =========
sobj_kera |>
  filter(orig.ident == 'infected') |>
  FeaturePlot(features = c('Trpv3','Ccl20'), blend = T, order = T,
              blend.threshold = 0, pt.size = 1)

Idents(sobj_kera) <- 'seurat_clusters'

### critial: define Trpv3-hi kera ----------------
sobj_kera <- sobj_kera |> mutate(trpv3_status = case_when(
  seurat_clusters == 5 ~ 'Trpv3-high',
  .default = 'Trpv3-low'
))

### FIG: trpv3-hi-lo umap ==========
sobj_kera |> DimPlot(group.by = 'trpv3_status') +
  theme_jpub

publish_pdf('micefig2/mice_v3h_v3l_umap.pdf', width = 60)

## investigate V3-hi --------------
### compare SA vs PBS in Trpv3-hi -----
trpv.h.infect.deg <- sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  FindMarkers(ident.1 = 'infected', group.by = 'orig.ident',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  filter(!str_detect(gene, '^Gm|Rik$'))

trpv.h.infect.deg |> write_csv('mission/FPP/mice_v3h_SA-PBS_deg.csv')

trpv.h.infect.deg |>
  filter(gene %in% key_cytokine)

### compare SA vs PBS in Trpv3-lo -----
v3l.sa.deg <- sobj_kera |>
  filter(trpv3_status != 'Trpv3-high') |>
  FindMarkers(ident.1 = 'infected', group.by = 'orig.ident',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  filter(!str_detect(gene, '^Gm|Rik$'))

v3l.sa.deg |> write_csv('mission/FPP/zww_sa_mice/results/mice_v3l_SA-PBS_deg.csv')

v3l.sa.deg |>
  filter(gene %in% key_cytokine)

### compare infected Trpv3hi vs Trpv3lo KC -----
sa_v3hvl_deg <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  FindMarkers(ident.1 = 'Trpv3-high', group.by = 'trpv3_status',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  select(c(gene, avg_log2FC, p_val_adj)) |>
  filter(!str_detect(gene, '^Gm|Rik$'))

sa_v3hvl_deg |>
  write_csv('mission/FPP/mice_SA_v3h-v3l_deg.csv')

sa_v3hvl_deg |>
  filter(gene %in% key_cytokine)

### compare ctrl Trpv3hi vs Trpv3lo KC -----
ct_v3hvl_deg <- sobj_kera |>
  filter(orig.ident != 'infected') |>
  FindMarkers(ident.1 = 'Trpv3-high', group.by = 'trpv3_status',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  select(c(gene, avg_log2FC, p_val_adj)) |>
  filter(!str_detect(gene, '^Gm|Rik$'))

ct_v3hvl_deg |>
  write_csv('mission/FPP/zww_sa_mice/results/mice_ctrl_v3h-v3l_deg.csv')

ct_v3hvl_deg |>
  filter(gene %in% key_cytokine)

### changes of Trpv3 upon infection --------
quant_v3_on_infection <- function(x){
  sobj_kera |>
  FindMarkers(ident.1 = 'infected',
              group.by = 'orig.ident',
              features = 'Trpv3',
              subset.ident = x,
              logfc.threshold = 0)
}

Idents(sobj_kera) <- 'trpv3_status'

trpv3_changes_upon_infec <- sobj_kera$trpv3_status |>
  unique() |>
  map(quant_v3_on_infection, .progress = T) |>
  set_names(unique(sobj_kera$trpv3_status)) |>
  list_rbind(names_to = 'cluster')

trpv3_changes_upon_infec

trpv3_changes_upon_infec |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(cluster, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = -1, linetype = 'dashed') +
  scale_fill_manual(values = c('blue','grey','red')) +
  labs(title = str_wrap('Change of Trpv3 expression upon infection', width = 30),
       fill = 'Significance') +
  theme_jpub

### FIG: TRPV3 SLE-HC logfc violin ----------
sa_v3h_meta <- sobj_kera |>
  as_tibble() |>
  select(.cell, orig.ident, trpv3_status)

sa_v3h_meta <- sobj_kera |> 
  get_abundance_sc_long('Trpv3') |>
  left_join(sa_v3h_meta)

sa_v3pval <- tibble(group1 = 'infected', group2 = 'PBS',
                 y.position = 3.2, label = c('***','*'),
                 trpv3_status = c('Trpv3-high','Trpv3-low'))

sa_v3h_meta |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(orig.ident, .abundance_RNA, fill = orig.ident)) +
  geom_violin() +
  stat_pvalue_manual(data = sa_v3pval, label = 'label',
                    label.size = 2, inherit.aes = F) +
  facet_wrap(~trpv3_status, ncol = 1) +
  theme_pubr() +
  scale_fill_manual(values = c('blue','red')) +
  expand_limits(y = 3.4) +
  labs(x = 'Group', y = 'Normalized expression',
       title = 'Trpv3 expression', fill = 'Group') +
  theme_jpub

publish_pdf('micefig2/mice_v3hl_trpv3_logfc_sa_violin.pdf')

### examine mva pathway ----------
sobj_kera |>
  mutate(group = str_c(trpv3_status, '_', orig.ident)) |>
  DotPlot(cholesterol_path, group.by = 'group') +
  RotatedAxis() +
  labs(x = 'gene', y = 'group')

ggsave('mice_SA_Trpv3hi+lo_cholesterol_pathway_bubbleplot.pdf',
       width = 9, height = 3, units = 'in')

#### FIG: logfc of mva pathway in v3h & v3l kera ===========
Idents(sobj_kera) <- 'trpv3_status'

t3h_mva <- sobj_kera |>
  FindMarkers(features = kegg_mva,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-high',
              logfc.threshold = 0) |>
  as_tibble(rownames = 'gene')

t3l_mva <- sobj_kera |>
  FindMarkers(features = kegg_mva,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-low',
              logfc.threshold = 0) |>
  as_tibble(rownames = 'gene')

v3kera_mva <- list(Trpv3_high = t3h_mva, Trpv3_low = t3l_mva) |>
  bind_rows(.id = 'celltype') |>
  left_join(chole_p_tib, join_by(gene == chole_p_str)) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > .3 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS'))

v3kera_mva |>
  ggplot(aes(cholesterol_path, avg_log2FC, fill = type)) +
  geom_col() +
  facet_wrap(~celltype, ncol = 1) +
  labs(x = 'Gene', fill = 'Significance') +
  scale_fill_manual(values = c('blue','grey','red')) +
  theme_jpub +
  RotatedAxis()

publish_pdf('micefig/mice_v3h_SA-logfc_mva_barplot.pdf', width = 70)

#### FIG: mva pbs-SA logfc in dotplot =========
v3kera_mva |>
  ggplot(aes(celltype, gene,
             color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  labs(x = 'Cell type', fill = 'Gene') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  RotatedAxis()

publish_pdf('micefig/mice_v3h_SA-logfc_mva_bubble.pdf', width = 50)

#### FIG: SA v3h vs v3l mva logfc dotplot =========
sa_hvl <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  FindMarkers(features = kegg_mva,
              ident.1 = 'Trpv3-high',
              group.by = 'trpv3_status',
              logfc.threshold = 0,
              min.pct = 0) |>
  as_tibble(rownames = 'gene')

sa_hvl |>
  left_join(mva_fct) |>
  ggplot(aes(y = fct_rev(ordered), x = 'Trpv3-hi keratinocytes',
             color = avg_log2FC, size = -log(p_val_adj))) +
  geom_point() +
  labs(y = 'Gene', x = NULL) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(breaks = c(0,100,200)) +
  theme_jpub

publish_pdf('micefig/mice_SA_Trpv3hi-lo-logfc_mva_bubbleplot.pdf')

### examine interested cytokines ------
sobj_kera |>
  mutate(group = str_c(trpv3_status, '_', orig.ident)) |>
  DotPlot(key_cytokine, group.by = 'group') +
  RotatedAxis() +
  labs(x = 'gene', y = 'group')

ggsave('mice_SA_v3h-v3l-pbs-SA_cytokine_bubbleplot.pdf',
       width = 6, height = 3, units = 'in')

#### FIG: logfc of cytokines in v3h/v3l kera upon SLE ==========
Idents(sobj_kera) <- 'trpv3_status'

t3h_key <- sobj_kera |>
  FindMarkers(features = key_cytokine,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-high',
              logfc.threshold = 0) |>
  as_tibble(rownames = 'gene')

t3l_key <- sobj_kera |>
  FindMarkers(features = key_cytokine,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-low',
              logfc.threshold = 0,min.pct = 0) |>
  as_tibble(rownames = 'gene')

list(Trpv3_high = t3h_key, Trpv3_low = t3l_key) |>
  bind_rows(.id = 'celltype') |>
  add_case(gene = 'Csf2', celltype = 'Trpv3_low', avg_log2FC = 0) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > .3 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(celltype, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 0) +
  facet_wrap(~gene, scales = 'free_y') +
  labs(fill = 'Significance', x = 'Cell type') +
  scale_fill_manual(values = c('blue','grey','red')) +
  theme_jpub +
  RotatedAxis()

publish_pdf('micefig/mice_v3h-v3l_SA-log2fc_cytokine_barplot.pdf', width = 60)

#### examine cytokine only in infection ----------
inf_cls_cytokine <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  SetIdent(value = 'seurat_clusters') |>
  FindAllMarkers(features = key_cytokine)

inf_cls_cytokine |>
  filter(p_val_adj < .05) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  theme_pubr() +
  labs_pubr()

publish_pdf('micefig/mice_SA_kera_cytokine_bubble.pdf')

# back to v3hi kera in all skin-----------
sobj <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

v3h_barc <- sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  pull(.cell)

v3l_barc <- sobj_kera |>
  filter(trpv3_status == 'Trpv3-low') |>
  pull(.cell)

sobj <- sobj |>
  mutate(bill_fine = case_when(
    .cell %in% v3h_barc ~ 'Trpv3-hi-Keratinocytes',
    .cell %in% v3l_barc ~ 'Trpv3-lo-Keratinocytes',
    manual_main == 'Keratinocytes' ~ 'bad-Kera',
    .default = manual_main
  ))

sobj %<>%
  mutate(bill_fine = case_when(
    str_detect(bill_fine, 'hi-K') ~ 'V3-hi KC',
    str_detect(bill_fine, 'lo-K') ~ 'V3-lo KC',
    .default = bill_fine
  ))

sobj |> write_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

sobj <- sobj |> filter(bill_fine != 'bad-Kera')

sa_fine_type <- sobj$bill_fine |>
  unique()

## MVA ------------
### FIG: v3hi kera in SA skin mva dotplot =========
sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine',dot.scale = 3) +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in SA-infected mice skin') +
  RotatedAxis()

sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine',dot.scale = 3,
          cols = 'RdYlBu') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in SA-infected mice skin') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_SA_skin_mva_bubble.pdf', width = 70)

### FIG: v3hi kera in PBS skin mva dotplot =========
sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine',dot.scale = 3) +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in control group mice skin') +
  RotatedAxis()

sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in PBS group mice skin') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_pbs_skin_mva_bubble.pdf', width = 70)

### FIG: mva enrichment in SA fine type ============
Idents(sobj) <- 'bill_fine'

fine_mva_enrich <- sobj |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers(features = kegg_mva,
                 logfc.threshold = 0,
                 min.pct = 0)

fine_mva_enrich |>
  mutate(cluster = fct_reorder(cluster, avg_log2FC),
         avg_log2FC = ifelse(p_val_adj == 1, NA, avg_log2FC)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  scale_size(breaks = c(0,100,200)) +
  #theme_jpub +
  labs(x = 'Gene', y = 'Cell type') +
  theme_pubr(legend = 'right') +
  RotatedAxis()

publish_pdf('micefig2/mice_SA_skin_mva_enrich.pdf', width = 70)

### FIG: mva change after SA infection in skin ==========
fine_mva_change <- sa_fine_type |>
  map(\(x)FindMarkers(sobj,
                      group.by = 'orig.ident',
                      ident.1 = 'infected',
                      subset.ident = x,
                      features = kegg_mva,
                      logfc.threshold = 0,
                      min.pct = 0) |> as_tibble(rownames = 'gene'),
      .progress = T) |>
  set_names(sa_fine_type) |>
  list_rbind(names_to = 'cluster')

g1 <- fine_mva_change |>
  mutate(cluster = as.character(cluster)) |>
  left_join(mva_fct) |>
  ggplot(aes(y = cluster, x = ordered,
             size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_pubr(legend = 'right') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Differential expression of MVA pathway: SA vs PBS') +
  RotatedAxis()

g1
g1 + theme_jpub +
  RotatedAxis() +
  scale_size(range = c(0,4))

publish_pdf('mission/FPP/micefig2/mice_skin_SA-PBS_mva_fc.pdf', width = 70)

## cytokine -----------
### FIG: cytokine expr in SA skin ===========
sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine') +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in SA-infected mice skin') +
  RotatedAxis()

sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in SA-infected mice skin') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_SA_skin_cytk_bubble.pdf', width = 60)

### FIG: cytokine expr in ctrl skin ===========
sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine') +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in control group mice skin') +
  RotatedAxis()

sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in PBS group mice skin') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_pbs_skin_cytk_bubble.pdf', width = 60)

### FIG: cytokine SA-PBS logfc in fine skin =========
sa_fine_fc_cytk <- sa_fine_type |>
  map(\(x)sobj |> FindMarkers(group.by = 'orig.ident',
                    ident.1 = 'infected',
                    features = key_cytokine,
                    logfc.threshold = 0,
                    min.pct = 0,
                    subset.ident = x) |> as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(sa_fine_type) |>
  list_rbind(names_to = 'cluster')

g1 <- sa_fine_fc_cytk |>
  mutate(cluster = as.character(cluster)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_pubr(legend = 'right') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Differential expression of cytokine: SA vs PBS') +
  RotatedAxis()

g1
g1 + theme_jpub +
  RotatedAxis() +
  scale_size(range = c(0,5))

publish_pdf('mission/FPP/micefig2/mice_skin_SA-PBS_cytk_fc.pdf', width = 60)

### FIG: cytokine enrich in SA skin ===========
fine_cytk_enrich <- sobj |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers(features = key_cytokine,
                 logfc.threshold = 0,
                 min.pct = 0)

fine_cytk_enrich |>
  mutate(cluster = as.character(cluster),
         avg_log2FC = ifelse(p_val_adj == 1, NA, avg_log2FC)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  scale_size(breaks = c(0,75,150)) +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type') +
  RotatedAxis()

publish_pdf('micefig2/mice_SA_skin_cytk_enrich.pdf', width = 60)

## save meta data for cell frac analysis ------
sobj_kera@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('mission/FPP/aureus_kera_meta.csv')

sobj@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('mission/FPP/aureus_skin_meta.csv')

# DC ---------
sobj |>
  DimPlot(group.by = 'manual_main')

sobj.dc <- sobj |>
  filter(bill_fine == 'DC')

sobj.dc %<>% quick_process_seurat(skip_norm = T, leiden = F)

sobj.dc |>
  DimPlot(group.by = 'manual_main')

dc.marker <- list(
  cDC1 = c('Itgae'),
  cDC2 = c('Itgam','Sirpa'),
  LC = c('Cd24a'),
  DN.DC = c('Xcr1'),
  MP = c('Fcgr1','Mertk')
)

sobj.dc |> DotPlot(dc.marker, cluster.idents = T)

sobj.dc |> FindAllMarkers(features = list_c(dc.marker), only.pos = T) |>
  filter(p_val_adj < .05)

immg <- celldex::ImmGenData()

immg.dc <- immg |>
  filter(label.main == 'DC')

sobj.dc %<>%
  mark_cell_type_singler(immg.dc, fine_label = T, new_label = 'immg.fine')

sobj.dc %<>%
  mutate(bill_fine = case_when(seurat_clusters == 5 ~ 'DN dDC',
                              seurat_clusters == 4 ~ 'Langerhans cells',
                                .default = 'cDC2'))

## DC subtype fraction changes ----------
sobj.dc |>
  as_tibble() |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  calc_frac_conf_on_grouped_count(group, bill_fine) |>
  ggplot(aes(group, fraction, color = group,
             ymax = conf.high, ymin = conf.low)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  facet_wrap(~bill_fine) +
  theme_pubr() +
  scale_color_hue(direction = -1) +
  labs(title = 'DC subsets proportion changes in skin SA infection')

sobj.dc |>
  as_tibble() |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  test_on_grouped_count(group, bill_fine)

sobj.dc |>
  mutate(subgroup = str_c(bill_fine, '_', orig.ident)) |>
  DotPlot(c('Ccr6','Mgl2'), group.by = 'subgroup') +
  labs(x = 'Gene', y = 'Cell type group')

sobj.dc |>
  filter(bill_fine == 'cDC2') |>
  VlnPlot('Mgl2', group.by = 'orig.ident', pt.size = 0)

sobj.dc |>
  filter(bill_fine == 'cDC2') |>
  FindMarkers(features = 'Mgl2', group.by = 'orig.ident', ident.1 = 'infected')

## CCR6 in cDC2? -------
immg.dc |>
  filter(.feature == 'Ccr6') |>
  ggplot(aes(label.fine, logcounts, color = label.fine)) +
  stat_summary() +
  stat_summary(geom = 'col', fill = 'white') +
  theme_pubr(legend = 'none') +
  coord_flip() +
  labs(title = 'Immgen: mouse DC expression of Ccr6')


# other immune cells -------
sobj |>
  DimPlot(group.by = 'manual_main', cols = 'Paired')

sobj |>
  as_tibble() |>
  filter(str_detect(manual_main, 'NK|T|Macro|DC|Gran')) |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  calc_frac_conf_on_grouped_count(group, manual_main) |>
  ggplot(aes(group, fraction, color = group,
             ymax = conf.high, ymin = conf.low)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  facet_wrap(~manual_main, scales = 'free_y') +
  theme_pubr() +
  scale_color_hue(direction = -1) +
  labs(title = 'Immune cell proportion changes in skin SA infection',
       y = 'Fraction in all immune cells')

## T cells ---------
sobj.ct <- sobj |>
  filter(manual_main %in% c('T cells', 'NK cells'))

sobj.ct %<>% quick_process_seurat(skip_norm = T)

sobj.ct %<>% mark_cell_type_singler(immg, fine_label = T, new_label = 'immg.fine')

mmur <- celldex::MouseRNAseqData()

sobj.ct %<>%
  mark_cell_type_singler(mmur, fine_label = T, new_label = 'mmur.fine')

sobj.ct %<>% filter(seurat_clusters != 7)

sobj.ct |> DotPlot(c('Cd3d','Il22','Gata3','Ccr4','Il17a','Il17f','Nkg7','Tnf','Ifng'))

sobj.ct |> VlnPlot(c('Ifng','Il22','Gata3','Il17a'), group.by = 'orig.ident',
                   ncol = 2)

sobj.ct |>
  FindMarkers(features = c('Il22','Gata3','Il17a','Ifng'),
              group.by = 'orig.ident', ident.1 = 'infected')

sobj.ct %<>%
  mutate(bill_fine = case_when(seurat_clusters == 4 ~ 'Th17 cells',
                               seurat_clusters %in% c(3,6) ~ 'Th2 cells',
                               .default = 'NK cells'))

sobj.ct |>
  as_tibble() |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  calc_frac_conf_on_grouped_count(group, bill_fine) |>
  ggplot(aes(group, fraction, color = group,
             ymax = conf.high, ymin = conf.low)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  facet_wrap(~bill_fine, scales = 'free_y') +
  theme_pubr() +
  scale_color_hue(direction = -1) +
  labs(title = 'NK/T cell proportion changes in skin SA infection',
       y = 'Fraction in NK/T cells')

# ATP receptor? -----
p2r.list <- rownames(sobj) |> str_subset('^P2')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(p2r.list, cluster.idents = T) +
  RotatedAxis() +
  ggtitle('ATP receptor in mice skin')

sobj_kera |>
  filter(orig.ident != 'PBS') |>
  DotPlot(p2r.list) +
  RotatedAxis()

sobj_kera |>
  filter(orig.ident != 'PBS') |>
  FindMarkers(features = p2r.list, group.by = 'trpv3_status', ident.1 = 'Trpv3-high')

# IFN receptor? -------
ifnr.list <- rownames(sobj) |> str_subset('^Ifn.r')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(c(ifnr.list, 'Il10rb')) +
  ggtitle('IFN receptor in mice skin')

ifnl.listm <- rownames(sobj) |> str_subset('Ifnl\\d')

sobj |>
  filter(orig.ident != 'PBS') |>
  DotPlot(ifnl.listm) +
  ggtitle('IFN lamdba expression in infected skin')

sobj |>
  mutate(subgroup = str_c(bill_fine, '_', orig.ident)) |>
  DotPlot(ifnl.listm, group.by = 'subgroup') +
  ggtitle('IFN lamdba expression in mice skin')

isg.listm <- rownames(sobj) |> str_subset('Ifitm')

isg.listm <- c('Oas2','Isg15')

sobj |>
  mutate(subgroup = str_c(orig.ident, '_', bill_fine)) |>
  DotPlot(isg.listm, group.by = 'subgroup') +
  ggtitle('ISG expression in infected skin')

# TLR? -------
tlr.listm <- rownames(sobj) |> str_subset('Tlr')

tlr.listm

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot('Sting1') +
  ggtitle('TLR expression in mice skin')

# RLR?
rlr.listh <- c('DDX58','IFIH1','DHX58')
rlr.listm <- rlr.listh |> str_to_title()

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(rlr.listm) +
  ggtitle('RLR expression in mice skin')

# inflammasome pathway -----------
nlrp3i.listm <- c('NLRP1A','NLRP1B','NLRP3','PYCARD','CASP1','IL1B','IL18',
                  'IL33','NLRP6','NLRP12','NLRC4','AIM2','MEFV','Ifi204') |>
  str_to_title()

sobj |>
  mutate(subgroup = str_c(bill_fine, '-', orig.ident)) |>
  DotPlot(nlrp3i.listm, group.by = 'subgroup') +
  ggtitle('Inflammasome pathway in SA infection') +
  RotatedAxis()

